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History-dependence and ageing are studied in the low-temperature glass phase of a long-range 
0^ , periodic Josephson array. This model is characterized by two parameters, the number of wires (2JV) 

■ and the flux per unit strip (a); in the limit N — *• oo and fixed a <C 1 the dynamics of the model are 

described by the set of coupled integral equations, which coincide with those for the p = 4 disordered 
spherical model. Below the glass transition we have solved these equations numerically in a number 
of different regimes. 

We observe power-law ageing after a fast quench with an exponent that decreases rapidly with 
temperature. After slow cooling to a not-too-low temperature, we see ageing characterized by the 
appearance of a new time scale which has a power law dependence on the cooling rate. By contrast if 
the array is cooled slowly to very low temperatures the ageing disappears. The physical consequences 
of these results in different cooling regimes are discussed for future experiment. 

We also study the structure of the phase space in the low temperature glassy regime. Analytically 
we expect an exponential number of metastable states just below the glass transition temperature 
with vanishing mutual overlap, and numerical results indicate that this scenario remains valid down 
to zero temperature. Thus in this array there is no further subdivision of metastable states. We 
also investigate the probability to evolve to different states given a starting overlap, and our results 
suggest a broad distribution of barriers. 
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I. INTRODUCTION 



a 

o. 

O , The essential features underlying the physics of glass formation remain elusive. A glassy system has a "memory" of 
its past history; it breaks ergodicity without thermodynamic selection of a unique state thereby defying description 
y—\ ' by the standard Gibbs methods. Though glass formation in the absence of disorder is a widespread phenomenon, 
. it remains poorly understood even at the mean-field level. Several infinite-range models for disorder-free glassiness 
have been proposed;EJ the majority were studied via a mapping to intrinsically random systems.! 2 ! This approach relies 
on the observation that certain infinite-range periodic and random models have identical leading order (in 1/N) 
high-temperature expansions. This correspondance is assumed to persist at low temperatures. Crystalline phase(s), 
which could appear in a periodic model, are missed by such a treatment; however the usual working assumption is 
that they are dynamically iiiaccessible. Recently a glass transition in a periodic, long-range Josephson array has been 
identified and characterizecatil using a generalized Thouless-Anderson-Palmer (TAP) approach.B This model can be 
studied directly (without a mapping to a disordered system) because it has two expansion parameters: the first is the 
usual l/N inherent to a long-range model; the second is a, the flux per unit strip, where the condition a <C 1 further 
restricts the number of contributing diagrams to be a small subset of all those appearing to O(jj-). The existence of 
a as an expansion parameter suggests that the crystalline phase(s) are inaccessible (or do not exist) for small values 
of a. We note that for a = 1 there exist non-trivial ground-states for an array with prime N which are certainly 
missed by perturbative methods. In this particular case, any sequence of complex numbers (representing the order 
parameters associated with each wire) of unit magnitude whose discrete Fourier transform has a flat spectrum yields 
a ground-state.tj This mathematical problem has been studied extensively in the context of number theory, beginning 
with the work of Gauss.D 

Another key feature of the periodic Josephson array is that it can be constructed in the laboratory; the physical 
requirements imposed by the theory on its realization have been discussed elsewhere.cl Measurement of the a.c. 
response ta_a small time- varying field in the fabricated array should probe its diverging relaxation time at the glass 
transition.Qa In this paper we study the behavior of this periodic network in its low-temperature non-ergodic regime. 
We find both history-dependence and ageing, and the physical consequences of these results for different cooling 
regimes are discussed for future experiment. 
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Fig 1. The phase diagram of the array (inset). Here Tq, T m and T5 are the temperatures associated with the dynamical 
instability, "superheating" , and a speculated equilibrium transition as discussed in the text. 



The system of our study is a stack of two mutually perpendicular sets of N parallel thin superconducting wires 
with Josephson junctions at each node (Figure 1) that is placed in an external tranverse field H . The classical 
thermodynamic variables of this array are the 2N superconducting phases associated with each wire. In the absence 
of an external field the phase differences would be zero at each junction, but this is not possible for finite H and the 
phases are thus frustrated. Here we assume that the Josephson couplings are sufficiently small so that the induced 
fields are negligible in comparison with H, an important point for the experimental realization of this network^ We 
can therefore describe the array by the Hamiltonian 

2N 

— ^ ^ 5 m j7mn'5n (1) 
m,n 

where J m n is the coupling matrix 

with Jj^ = -jf= exp(2iriajk / N) and 1 < (j, k)<N where j(k) is the index of the horizontal (vertical) wires; s m = e l ^ m 
where the <f> m are the superconducting phases of the 2N wires. Here we have introduced the flux per unit strip, 
a = NHl 2 /$ , where I is the inter-node spacing and $0 is the A ux quantum; the normalization has been chosen so 
that Tq does not scale with N. 

Because every horizontal (vertical) wire is linked to eyery vertical (horizontal) wire, the connectivity in this model is 
high (N) and it is accessible to a mean-field treatment For a < 1 there exist an extensive number of metastable 
solutions separated by barriers that scaletl with N. A similar long-range network with disorder was previously found 
to display a spin glass transitionB for a 3> jj] in the absence of short-range phase coherence between wires (a ^> 1), 
it was mapped onto the Sherrington-Kirkpatrick modelJl3 Physically this glassy behavior occurs because the phase 
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differences associated with the couplings, Jjk, acquire random values and fill the interval (0,2n) uniformly. For the 
periodic case, this condition is satisfied in the "incommensurate window" <§; a ^ 1; here the magnetic unit cell 
is larger than the system size so that the simple "crystalline" phase is inaccessible B There are thus no special field 
values where the number of low-temperature solutions are not extensive, in contrast to the situation for a > 1 . Of 
course, we cannot exclude the possibility that at low temperatures a single low-lying state appears; however it seems 
very likely that even in this case it will be irrelevant because the system is trapped in one of the many metastablc 
glassy states. Furthermore we have not observed such states in our direct numerical simulations of a = 1 array. 

In the thermodynamic Jimit of N — ► oo (with fixed area), the system is,-accessible to a static and dynamical 
high-temperature study;Ba this was performed using a modified TAP method.El This will be discussed in more detail 
when the coupled dynamical equations are derived in Section II; at temperatures above the dynamic instability the 
system is still in equilibrium and these equations can be solved analytically. The glass transition is characterized by 
a diverging relaxation time and a jump in the Edwards-Anderson order parameter, and does not coincide with an 
accompanying static transition. This calculated glass transition temperature (Tq) assumes that cooling from high 
temperatures occurs infinitesmally slowly; in practice, the system's effective glass transition temperature (T G > Tq) 
will correspond to the system's failure to equilibrate and will be a function of its finite cooling time-scale. For the sake 
of notational simplicity, we will not distinguish between the effective and the "adiabatic" glass transition temperature 
in the text that follows unless explicitly stated. 

Below the glass transition temperature (Tq), the system can no longer equilibrate. In contrast to the situation 
for T > Tg, the response and the correlation functions are not a function of time differences but rather are history- 
dependent; furthermore we can no longer relate them with the fluctuation-dissipation theorem. This makes the 
solution of these equations rather complicated. The low-temperature state of the system is determined by the sample 
history; here we shall consider the simplest scenario in which the temperature is reduced linearly from an initial T 
slightly above Tq to a final temperature T/(< Tq) in a time t c . The final state of the system at time t > t c is 
controlled by both time-scales t and t c , and we expect qualitatively different behavior depending on their relative 
magnitude. Here we shall focus on the limiting cases t 3> t c and t = t c . 

Physically in the limit of very fast cooling (t 3> t c ), we expect the state immediately following the quench to have 
minimal overlap with that observed at time t so that D U '^q = 0. Assuming that the solutions of the dynamical 
equations are well-defined in the limit when the microscopic time-scale goes to zero, we have only two time-scales in 
the problem: t and t' . In this case D tt i may depend only on their ratio t' jt on dimensional grounds: D t f — d(K). 
Since d(0) = we expect lim t // t _ i0 d(j) = (y) 7 Z)o, and this functional form of D u > remainstiHl3 at larger values of 
(£). Similarly in this fast-cooling regime we expect lim t7t G«, = 

F (t) 7 G o> so that in thc limit of small (*7*) 
there is a quasi-FDT relation (G cx D'), which remains vahcfiEl over the full range < (t'/t) < 1. 

The other limiting case is that of slow cooling (t = t c ), where we expect that the correlation and response functions 

will have their only time-dependence through the specified cooling process r(t) where r = ^ T 9 r ~ T ^J : D t v = d[r(t), r(t')] 

and Gu' — dT }^, ^ g[r(t), r(t')]. If this conjecture is correct, it implies that the system has a non-decaying memory 
of its past history. Of course this assumes again that the solutions to the dynamical .equations are well-defined in 
the limit of vanishing microscopic time-scales. For the Sherrington-Kirkpatrick model,EJ this "adiabatic ansatz" for 
linear r(t) leadstS to the solutions d(t,t') = t# and g(tjQ = 2r t /, which result in a field-cooled susceptibility in 
qualitative agreement with a previous equilibrium result x3 However application of the same method to the p = 4 
Ising model leads to a delta-function solution for the response function;Ej this form is worrisome given the original 
adiabatic nature of the ansatz. However, in principle the response could be, broadened to be long compared to the 
microscopic time-scale but short compared to the final measurement time.Ej When we started our numerical study 
of this problem, our expectations were therefore that the only time scale in the slow cooling process is set by the 
cooling rate, t* = (t|) _1 5 and that the response and the correlation functions depend only on the ratios t/t* and t/t* 
as indicated in the t = t c column of Table 1. This hypothesis turned out to be only partially correct. Instead we 
have found that indeed just below Tq the correlation and response functions depend only on the ratios t/t* and t/t* 

but t* is a new time-scale t* = (jrf) 1 ''^ that emerges in this cooling regime; however there is evidence that thc 
exponent n vanishes at low temperatures, indicating that "adiabatic cooling" is recovered. We note that a finite value 
of rj implies a decaying memory and hence ageing, and thus the disappearance of n means that the system recovers 
its memory and thus "rejuvenates" at low temperatures. 

Our expectations for the three cooling regimes are summarized in Table I. In Section II. we derive the dynamical 
equations for the array in the low-field (a <C 1) regime. The system goes out of equilibrium at a temperature Tq 
(determined by the cooling rate), and we study its low-temperature (T < Tq) behavior numerically in all three cooling 
regimes; our results are presented in Section III. 

We also probe the structure of phase space in this glassy regime (Section HID) where the periodic array is stuck in one 
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of an extensive number of metastable states.0 What happens to these states when the temperature is decreased? Each 
state, formed at To, may evolve smoothly; alternatively it can further subdivide as in the case of the (disordered) 
Sherrington-Kirkpatrick model. We find that in this periodic array there is no further subdivision down to zero 
temperatures, and that the number of metastable states remains exponentially large. We also learn that though these 
states are equidistant in phase space, their associated basins of attraction are rather complicated. We conclude with 
a discussion of the physical consequences for experiment (Section IV), a brief summary and with projects for future 
work. 

TABLE I. Three distinct regimes and their associated response (D tt i and correlation (G U ') functions when the system is 
quenched out-of-equilibrium (t < tit); here t, t r and t c are the measurement, relaxation and cooling times respectively, and t* 
is a time-scale which is a function of the cooling rate (4? ) • 

Du> = (?) ' D m At'=d(£,£) 

G u ,=xD' tt , ??? G«, = £ fl (£,£) 

**=/(£) 
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II. THE DYNAMICAL EQUATIONS 



The time-dependent properties of a glass are particularly distinctive, and thus we begin our study with a discussion 
of the self-consistent equations describing the dynamics of the periodic Josephson array. For the sake of convenience, 
we shall refer to the network model described by the Hamiltonian (|) in "magnetic" language so that each "spin" 
s m = e*^ m is associated with the phase of a superconducting wire labelled by "site index" m. We presume that the 
system's long-time behavior can be described by that of soft spins with pure relaxational modes, since intrinsic single- 
spin dynamics should be irrelevant on these time-scales. More specifically we introduce a "constraining potential" 
V{si) — Vb ( | s | 2 — l) 2 at each site that restricts the magnitude of each spin, |sj| » 1, and then assume the equations 
of motion 

lcHH + V) 

nsi = -y — 7^ — + Q (3) 

where t& is a microscopic time-scale; the effects of a coupled heat bath are represented by a time-varying Gaussian 
random field Q(t) with zero mean and variance 

(C i (t)C j (f)) = 2 n S ij S(t-l/) (4) 

where the brackets refer to an average over the noise. This Langevin equation (||) reproduces the dynamics of the 
overdamped Josephson array! with effective resistance R if r\, — ^3. ^ and Vq — > oo where Rq = ; here R — 
where Rq is the individual junction resistance since the circuit elements in the long-range array operate in parallel.H 
The physical quantities of interest here are the dynamical cumulants of the "spin" variables averaged over the noise. 
In particular, we would like to determine the time-dependant two-spin correlation function 

D mn (t,t') = (s m (t)s* n (t')) (5) 

and the linear response 

where again the angular brackets refer to an average over Q(t). We note that G mn (t,t') is a susceptibility with 
respect to a "theorist's" field h n conjugate to a spin s n = e** 1 ", and thus cannot be measured directly; however once 
known, G mn (t,t') can be.-used to determine the ac-response to a time-varying physical external field H(t) which is 
experimentally accessible.u 

Exploiting the fact that the Langevin noise Ci(t) is Gaussian and uncorrected in space and time, we can perform 
the noise peerages in (|J) and (^) using a functional integral formulation;E3EZl details of this technique can be found 
elsewhere.liS For example, the two-spin dynamical correlation function is given by the expression 

D mn {t,t') = (s m (t)s* n (t')) = J s m (t)s* n (t)expA[s,s]VsVs. (7) 

with the action 



A[s. s] = j ,11 s I T b S + L ^2LtZl \ + Tb § 2 + h.C. 



(8) 



where we have neglected the term arising from the functional Jacobian since it does not affect the system's long-time 
behavior provided that causality is maintained^ the latter is ensured by the diagrammatic treatment described below. 
The auxiliary field s in a correlation functionlij acts like a response term so that G mn (t,t') = (s m (t)s n (t')) can 

be determined by an expression analogous to (Q). 
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(d) 



-e- 











Fig 2. The diagrammatic expansion (to leading order in 1/N) for the (a) response and the (b) correlation functions. 

Here the coupling matrix J is represented by dashed lines, V is the single-site four-spin vertex shown as a diamond, Gu is the 

single-site irreducible response (correlation) function, G (D) is the full response (correlation) function, and Gj, (-D(,) is the single-site bare 
response (correlation) function. The leading self-energy contributions to G u and D u (in a -C 1) are displayed in (c) and (d) respectively. 



More specifically, we calculate the time-dependent correlation and response functions, D mn (t,t') and G mn (t,t') 
respectively, by resumming the terms in the ^ expansion of the generating functional associated with (||) which are 
leading order in ; this is a dynamical generalizationtrEl of the Thouless- Anderson-Palmer approachB where feedback 
through single-site Onsager terms play a crucial role. It is convenient to determine these response and correlation 
functions initially in the frequency-domain. The leading diagrams (in for G w are shown in Figure 2a; in the 
thermodynamic limit (limiV — > oo with fixed array area) the subleading corrections to the response can be safely 
neglectedfl Summing the geometric series shown in Figure 2a, we obtain 

G w = - =- (9) 

G^-/? 2 (JJt)G w 

for the response function connecting wires of the same type (horizontal/ vertical); here we consider the local Green's 
functions G w that is irreducible with respect to the Ju lines. The matrix (JJ^)y depends only on the "distance" 
i — j and acquires a simple form in Fourier space lJJ')p = ( J§ /a)9(air — \p\); therefore in this representation the 
low-frequency contribution to the response functionu is 

(0) _ f dp 6(aTr-\p\) 9(\p\-air) 

J (2ir) G -^(mi Guj + GZ 1 1 J 

which is of primary interest in this problem. We note that the static limit (uj — 0) of TG" 1 coincides with the locator, 
A(T), discussed in previous work on this array ;i in the absence of feedback terms Go =1. We also see in ( |l0| ) 
that there is a static instability at Gq 1 = G c = ^= in the vicinity of To = 77=- For 9 = < " T t^°' > a ^ feedback 

terms are negligible and Gj 1 = Gj, (lu) = f3A{T) — iujTb where t\, is a microscopic time-scale. However at lower 
temperatures, < -^/a, feedback effects introduce an additional frequency-dependent part of the local response 
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G w = G h {u>) + S w (11) 

as a self-energy (E u such that £o = 0) since we have chosen our normalization so that Gq 1 = (3A(T). The leading 
contribution in a to the self-energy E w is displayed in Figure 2b. We note that the self-energy contains a term 
6 = /3oJ^GJ where G is the full response function and Pq = -jr, and it is therefore convenient to write an equation 

for G. We shall use the identity 

G = (G-,SgjtGJ) + G, G=/?2jtGJ (12) 

and note that, since the matrix JJ^ has only two eigenvalues, and Tq , the expression (G— (3q J^GJ) = (1— 0$ jtj)G = 
(1— /3gJ^J)G projects only the zero eigenvalue contribution of JJ^ in the denominator of (^). We note that G is 
diagonal so that the order of matrix multiplication is not important. G is equivalent to the long-time part of the 
response function and is proportional to a (cf. (|l0|). Higher order terms in the self-energy will thus contain higher 
powers of a, and thus we conjecture that for a <C 1 these terms can be neglected; this will be justified later by the 
structure of the resulting solutions for G and D. 

The retardation of the self-energy terms (cf. Fig. 2b and 2d) significantly affectpthe long-term behavior of the system. 
This feedback results in a diverging relaxation time and a dynamical instability-! at a temperature Tq = (1 + Qq)To 

where Qq = — Ir^-a 1 / 4 . 

For temperatures T > Tq the response and correlation functions depend only on the time-differences, and thus can 
be related by the Fluctuation-Dissipation Theorem. However for T < Tq this relation no longer holds and we must 
determine the correlation function independently. Analogous to our treatment above, we introduce D and D as in 
equation ([l2|): 

D = (D-/3gj+DJ) + D, D =$gjtDJ. (13) 

where again we use (D— [3q J^DJ) = (1— /3qJTJ)D. The leading diagrams (in -k) for D u are shown in Figure 2c; the 
low-frequency contribution to this correlation function is 

J {2k) G^ - i^ Gui " G-J ^G u [ ) 

where n is an on-site self-energy contribution to D u . For long times and (|<)a<Cl 

a 

Gti,t 2 = — — —q tt; (15) 



and 



U ti,t, —ITTi ^ (I 6 ) 



where a = f3(A — ^) « 20 for 9 yfa and a k -| for — ^> *Ja, and A is the locator .□ Possible on-site self- 
energy terms, S and n, displayed to leading order in a in Figures 2b and 2d, simplify in the limit of a <C 1. In this 
regime, which we consider here, we note from ( p"5| ) and ( |l6| ) that the slowly decaying parts of G and D scale with a. 
Therefore the dominant self-energy contributions contain the minimal number of these functions (cf. Figs. 2c and 2d 
respectively). We can therefore consider only the leading order feedback terms 

3 * 2 

Eti,t2 = -^D tlt2 G t i,t2 (17) 

and 

n*M2 =4 £ U + 2 nSa,t2 (18) 



where G and D are as in ( |15|) and (|16|) respectively, and we have approximated the four-spin vertex by its static value 

r = -i. 

We now have expressions for the correlation and the response functions 
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a tl + 2r b ^ - 2St 1)ta ) G tlM = aS{h - t 2 ) (19) 



a tl +2r 6 J- -2E tl , t2 ) D tuta = / Dl M G t2iV dt' + 2r h G tlM . (20) 



dti * iLJ t 1 ,t2 J ^ti,*2 J -^ti.ta 

In order to solve these equations we also require a boundary condition for D(t,t), which follows from the definition 
( |l3| ) and the condition that Di,i(t, t) = 1 so that 

(l-a)D(t,t)+D(t,t) = l. (21) 

We only need to consider equation ([HJ) to order 0(a 1 ^ 4 ) which allows us to neglect self-energy corrections to D (cf. 
Figure 2d) that scale as a 3/4 in this case D{t, t) w D {0 \t,t) then i)(t,t) w 1 - D^(t,t). In the vicinity of the glass 
transition studied here, — ~ a 1 / 4 3> a/cJ, we have@ D$ = G^°2 w (1 + 0) and therefore 

Dt,t = -©f (22) 



This boundary condition and equations (|19| ) and (|20|) form a closed system of equations. 
We can simplify ( |l9| ) and ( |20| ) by defining the parameters 

t = 2na- s,i = -0oT 1 / 4 a tl = a tl a 3/4 tl (23) 

and the rescaled functions 

t) tu t,= ^%J^-,^f] (24) 
\t v t t v t2 l 



1/4 1/^1 ^2 



to obtain the equations 

(«ti + Ml) ff(*i.*a) ~ 303 J* <P(t 1 ,t)g(t 1 JW,t 2 )0t'dt!= 5{h - t 2 ) (26) 

(ot, + J r ) dfa.ia) - 302 /J 1 (Pih^gih^git'^dldf 

-9lf* 2 d 3 (t 1 ,t')g(t 2 ,t')e t ,dt'-2g(t 1 ,t 2 ) = (27) 

where we have used the leading order (in a -C 1) self-energy terms, ( |l7| ) and (|l8|). The boundary conditions, d(t, t) = 1 
and 
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7(iii*2)| tl _ >t + =1 ^-d(h,t 2 ] 



(28) 



yields another relation 

d t = l + 40 3 / d 3 {t,t')g(t,t')6' t dt' (29) 



which completes the closed set of dynamical equations for the system. From (pa) we see that the response and the 
correlation functions are proportional to a 1 / 4 , thus justifying our earlier assumption that higher-order terms in the 
self-energies are smaller in a. 

We pete that (|26|) , (p7|) and (^9|) are identical to the self-consistent equations for the p = 4 disordered spherical 
model,E3 which is particularly interesting since there is no intrinsic disorder in the periodic array described by (^) in 
this study. For T > Tq, the dynamical equations, (p(f), (|27j ) and (|29|), can be collapsed into one equation using the 
fluctuation-dissipation theorem and time-translation invariance; as expected from known results on p > 2 spherical 
models,o there is a dynamical instability that is unaccompanied by a static transition!] The introduction of higher 
order (in a) self-energy terms leads to dynamical equations for p disordered spherical models with several coexisting 
values of p (i.e. p = 4,6.. where p > 4 and p is evpia); thus this periodic array is the realization of a "simple spin 
glass" that was previously a theorist's model systemcil without a clear experimental counterpart. We note again that 
the disordered array previously studiedo is qualitatively similar to the p — 2 disordered XY model, and it is intriguing 
that the absence/presence of disorder in the array leads to the absence/presence of coinciding static and dynamical 
transitions. 
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III. THE COOLING REGIMES 



The self-consistent dynamical equations, (|26|), ( |27| ) and (|29|), can be simplified substantially^ for temperatures 
T > Tq where the response and the correlation functions depend only on the time-differences and can be related by 
the Fluctuation-Dissipation Theorem 

g{h - 1 2 ) = - arf( d~ t2 W i - 1 2 ). (30) 

Here mode-coupling theory pis exact because of the absence of a length-scale in the interaction; following previous 
discussions of this approached we obtainedo the closed form equation 

1 + r r^X <d , ( 31 ) 

o 1 iuj iuj Jq e IUJ td^(t)9^dt 

which results in a scaling form for the relaxation time 

t R = t T- u (32) 

where r = 9< g~ 9 , v = 1.765 and t is a constant that can be determined numerically .§ This identification of the 
glass transition at 9q assumes infinitesmally slow cooling from high temperatures; then the system will only fail to 
equilibrate when the relaxation rate diverges. 

In practice the periodic array will be subject to cooling at a finite-rate. In this case, the system will fall out of 
equilibrium at a reduced temperature {9 e J^ < 9q) when its time-scales associated with cooling and relaxation are 
equivalent, namely when the condition 

t T~ u (33) 



dr/dt 

is satisfied. (33) leads to the simple relation r 1+l/ = t Q (^) t a (^f), so that for a given finite cooling rate 4| the 



array will be "out of equilibrium" for times t > t B Q where 



4 s = t T-» = t \t ( ^ ) I- (34) 



d6 
Jtt 

and undergo an effective glass transition at 6c e ^ = ) wn i cn i s a function of its cooling rate. In the interest of 

simplicity, we will no longer distinguish between the system's effective and "adiabatic" dynamical instability; however 
here we wish to emphasize that their distinction is related to the finite nature of the cooling rate. 

We would now like to characterize the array's response and correlation functions at low temperatures (9 > 9 G ) in 
its glassy regime with the hope of identifying memory and history-dependent features. Unfortunately, in the absence 
of time-translational invariance and the fluctuation-dissipation theorem, the solution of the self-consistent dynamical 
equations derived in Section II. is difficult. As already described in the Introduction, the expected behavior of Gu' 
and Du> is qualitatively different for the limiting cases of fast (t ^S> t c ) and slow cooling (t — t c ) regimes; here the final 
temperature T/ < Tg is reached at time t c . Furthermore at the onset of this study, there were no clues about the 
physically interesting (and experimentally accessible) crossover region where t>t c . In order to clarify this situation, 
we now describe our numerical studies of the dynamical equations (|26|), ( p7| ) and ( |29j ) in the three distinct cooling 

o3/4 

regimes. Here, in the interest of simplifying our notation, we measure all reduced temperatures in units of 9q — 
so that 

9 -> 9 G 9. (35) 



A. The Fast Cooling Regime (t > t c ) 

In the limit of fast cooling all measurements are performed at times long after the system has reached its final 
temperature. Recently there has been p rocra ss in the characterization of infinite-range (disordered) spherical models 
with p-spin interactions in this regime pM3 these results for the p ~ 4 case are relevant for the periodic array. 
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Furthermore the predicted scaling forms of the response and the correlation functions lead to physical consequences 
that could be experimentally observable; this will be discussed in Section IV. After a fast quench, the correlation 
function is predictedy'fca to take the form 



t'-tp 
t-t 



d 



(36) 



where to = t c + t w is the time when the system reaches its final temperature (T < Tq] here t c and t w are the cooling 
and waiting times respectively. In Figure 3 we display du> plotted as a function of i = * for different cooling and 
measurement times where Of = 1.1 and < t < 0.96 corresponds to the "fast quench" regime. Short-time dynamics 
are observed in the range t > 0.98, and t < corresponds to the "time range" when the temperature was changing 
at a finite rate; this range vanishes in the limit of very fast cooling. In Figure 3, we see numerical support for the 
scaling of the correlation function implied by expression ([56]); we see there that with increasing t, the correlation 
functions corresponding to finite-cooling rate converge to that for infinitely fast cooling. A log-log plot of the dtt> 
for the "fastest" response (with the least "tail" for t < 0) is linear, and yields 7 = 0.43-ajad do = 0.8 for Of = 1.1, 
consistent with the numerical trends previously reported for the p-spm spherical models .l n H 12 l In the case of infinitely 
fast cooling to = t w + t c , where t w is the waiting time at the initial (high) temperature. Furthermore we found 
that the correlation function for finite-rate cooling is well-described by the scaling form ( [36] ) even for relatively slow 
cooling rates if to is adjusted so that it corresponds to the time at which the system falls out of equilibrium; now 
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Fig 3. Numerical evidence for scaling of the correlation function D tt i after a fast quench to reduced temperature 
9f = 1.1. The solid line corresponds to the case of an almost infinitely fast quench. The other curves correspond to a finite-cooling rate of 
^2 = 0.05 to the same reduced temperature with different measurement times t=200,300 and 400 (top to bottom); we observe that these 
finite-rate correlations converge to the fast quench case with increasing t. 



A generalized fluctuation-dissipation theorem (gFDT) 

9tt 



(37) 



has also been identifiedSS in this cooling limit for the range < t < 1. In Figure 4 the numerical results strongly 
support this proposal (note that t > 0.8 and t *~ correspond to short-time dynamics and to time-interval where 
the temperature is var.yj.ng); furthermore the value for x — 0.55 is in close agreement with that expected from the 



analytical treatmentErE3 {x = (p — 2)(1 — q) with p = 4 and q — |). We also note that the scaling form (36) appears 
to be robust to decreasing final temperature, though the specific values of do and 7 change (for Of — 1.3, do = 0.95 
and 7 = 0.2 in contrast to the values d = 0.8 and 7 = 0.43 for Of = 1.1). 
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(t- - t )/(t - g 

Fig 4. Numerical support for the generalized FDT relation after a fast quench. Here we plot x(t') = J 1 /^'*/^ as a function of 
t = (t' — to)/ (t — to) which is constant but less than unity for t < 0.8 where to corresponds to the time when there is no longer temperature 
variation; t > 0.8 corresponds to short-time dynamics. 



B. The Slow Cooling Regime (t = t c ) 

The other limiting case is that of slow cooling where the measurement is taken at the final cooling time (t = t c ). 
We expect t* , a time-scale that is a function of the cooling rate T = (ir), to play a key role in the time-dependence 

of the response and correlation functions; more specifically the latter will be a function of p- and ^ as displayed 
in Table I. where t* = f (4). The form of the functional dependence of t* on the cooling rate must be determined 
numerically. In Figure 5, we plot and gu< vs. Tt' where V — for Of = 1.1. We remark that in the absence 
of anomalous scaling (rj = 0), dtt' an d gw are simply functions of the reduced temperatures Ot and 0' t so that in this 
regime they can be plotted in terms of r = 6' t — 1. We observe that the characteristic width decreases with increasing 
cooling times suggesting the need for another functional form for T. In Figure 6, we present another scaling plot for 
the same data as in Figure 5 (Of = 1.1); this time we plot and gw vs. ft' where f = r 1- ' 7 with r\ = 0.19; now the 
curves fall on top of one another nicely. In this case du' and gtt> are not just functions of the reduced temperatures. 
Scaling plots are also presented for data associated with Of = 1.2 and Of — 1.3 in Figures 7 and 8; we infer from these 
results that 



(38) 

where rj decreases with decreasing temperature. More specifically we find that 77 = 0.19, 0.12 and 0.0 for Of = 1.1, 1.2 
and 1.3, and the associated scaling plots are displayed in Figures 6-8. In Figures 6 and 7 n 7^ and there is a 
new time-scale, t* , associated with the decay of the correlation functions. The scaling plots in Figure 8 indicate 
that r7 = at lower temperatures, suggesting that the array "rejuvenates" and improves its memory there. We 
note that 77 = in— the slow cooling regime of the Sherrington-Kirkpatrick modelji3 it is somewhat amusing that 
such ageing effectspJ'LJ well-known in spin glasses, should appear in a model for a periodic array. The observable 
consequences of this "forgetfulness" and subsequent "rejuvenation" will be discussed in Section III. Finally we note 
that at temperatures such that 77 = dd g t f ^ inconstant for a wide temperature range ressembling the slow cooling 
solution of the Sherrington-Kirkpatrick modelta. We also recall that there were no qualitative changes in the fast- 
cooling regime as a function of increasing Of (decreasing temperature). 
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Fig 5. Absence of the adiabatic limit in the slow cooling regime. Here we plot jr ^'gt'* ant ^ rS(*>' ) vs - r = — *g) = 

#(t') — 1 for different cooling rates (T = 0.3/40, 0.3/80 and 0.3/160); the characteristic width of the curves decreases with increasing cooling 
times. All dashed curves are d'(t,t') and the solid ones are g(t,t'), and the cooling time increases from bottom to top. 
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Fig 6. A scaling plot of i dd £ t 'f' ) vs r(f) = f(t' - t) for three different cooling rates (T = 0.3/40,0.3/80 and 0.3/160) where 
f = r 1— 11 with rj = 0.19; here 8f = 1.1 All dashed curves are d'(t, t') and the solid ones are g(t, t'). 
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Fig 7. A scaling plot of I ^gf 1 vs r(t') 



f(t' - t) for three different cooling rates (T = 0.3/40,0.3/80 and 0.3/160) where 



f = r 1 v with rj = 0.12; here 8f = 1.2. All dashed curves are d'(t, t') and the solid ones are g(t, t'). 
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Fig 8. A scaling plot of I dd ^P vs r(t') = r(t' - t G ) = 9{t') - 1 for three different cooling rates (r = 0.3/40, 0.3/80 and 0.3/160) 
indicating the appearance of the adiabatic limit at 9f = 1.3. All dashed curves are d'(t,t') and the solid ones are g(t,t'), and the cooling 
time increases from bottom to top. 
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C. The Crossover Cooling Regime (t > i c ) 



As indicated in Table I, we had no initial expectations for the form of d U ' and gtr in the intermediate cooling regime 
t ~> t c . In Figure 9 we display a number of correlation functions that display crossover behavior between the two 
limits of the fast and slow cooling cases. Here we see that at t > t c there exist at least two distinct regimes: t' < to 
and t' > to. In the first regime is approximately constant, whereas in the second case the correlation function 
varies more like a power-law ressembling the fast cooling case. In order to emphasize this point, we have included a 

scaling curve D{t') — do ffEf) where i = \zl° \ here x — 0.13 is parameter that effectively rescales to, and we have 
taken the parameters do — 0.8 and 7 = 0.43 from the fast-cooling case (see Figure 3) where 9f = 1.1. We note that 
the scaling curve and the correlation function (t — 432) agree very well for t > 0.3; we also see an intermediate regime 
< t < 0.3 where there is still curvature in D(t') but no agreement with the scaling form (^). The generalized FDT 
ratio x, displayed in Figure 10, is as previously observed for fast cooling (cf. Figure 4) for t — * ~*° > 0, but is no 

longer constant for t < (slow cooling). By contrast in Figure 11 we observe that dd g t ', t ^ is almost constant as in the 
slow-cooling regime discussed above (see Figs 5-7). This overall behavior suggests that for "crossover measurement 
times" (t > t c ), the system's behavior is a combination of that observed for both fast and slow cooling. As discussed 
in Section IIB, the behavior of the slow-cooling correlation function changes qualitatively when the final temperature 
is sufficiently low {Of > 1.3). This change also affects the correlation function in the crossover regime. As we observe 
in Figure 12 the ^ddt^ acquires a peak and, more importantly, it ceases to depend on the cooling rate T indicating 
that the system preserves its memory of the cooling process even as t — > 00 (with the ratio — maintained constant). 
We note that the plateau displayed in Figure 12 corresponds to that of d' shown in Figures 7 and 8 (in the slow 
cooling regime). The range t — t Zl° > corresponds to the fast cooling regime where ~ ^ ; the observed peak 
is needed to match the fast- and slow- cooling regimes. 




(f - t )/(t - t ) 



Fig 9. The correlation function d(t,t') in the crossover between fast and slow cooling regime for cooling rate V = 0.0125 

and t = 180, 232, 332, 432 and 6f = 1.1 (top to bottom). Note the presence of two different regimes evident for t ~ t c : (1) ~~g|r ~ 

~ ~ ~ t' ' —t 

constant for t < and (2) power-law variations at t > ressembling the fast cooling regime where t = . Note that at t becomes 

r — 1 

larger than t c the linear slow-cooling regime rapidly vanishes and that the correlation function converges to the scaling form (Bq) (dash-dot 
curve) D = d (jEf ) 7 with d = 0.8, 7 = 0.43 and x = -0.13. 
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Fig 10. The generalized FDT ratio x in the crossover regime between slow and fast cooling for T = 0.125, t = 432 
and 9f = 1.1. Note that in this plot t was sufficiently larger than t c and the modified FDT relation is maintained for most of the range 
(((' — to)/(t — to) > 0); furthermore x is as previously observed for fast cooling, but for (t' — *o)/(t — to) < the ratio is no longer constant. 
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Fig 11. The derivative of the correlation function, ^jP- for F = 0.004, t = 180 and 9 f = 1.1 in the crossover between 
slow and fast cooling. In this figure t ~ t c and we see the appearance of the behavior characteristic of slow cooling ( g'f ~ constant). 



15 



0.5 



1 1 

- e i = 


I 

0.85, 6, = 1.30 


1111,1111 

t c = 100, t = 200 . 


- t. = 




1 

I 


t c = 144, t = 288 . 


- 

- // 

/ / 






j 

1 1 


I I 1 I I I I 


i i i i 1 i i i i 



-1 -0.5 0.5 1 

(f - t )/(t - t ) 



Fig 12. The derivative of the correlation function, ^pg^ - for two cooling rates, V = 0.45/100, .45/144 with the same ratio 
= 2 and 9f = 1.3. We see that d ^g^ ^ is constant even in the slow cooling regime indicating "perfect" memory. 



D. Non-Monotonic Cooling and Structure of the Phase Space 



It is known that the structure of the metastable states in the (disordered) Sherrington-Kirkpatrick (SK) model is 
rather complex; the states are arranged on an ultrametric tree with many branch points. Here we provide evidence 
that the situation is rather different for the periodic array discussed here. In the glassy phase of the SK model, each 
metastable state subdivides and "daughter" states appear as the temperature is lowered. In this case application of 
a heat pulse removes some of this subdivision and results in a partial loss of memory. For the periodic array such 
behavior is not observed but, by contrast, a small heat pulse does not lead to any aftereffect. Of course a large heat 
pulse brings the system to above its effective glass transition temperature and thus erases all of its memory. In Figure 
13 we display the correlation functions for varying amplitudes of applied heat pulses. Here the system was cooled 
fast to low temperature (9 = 1.4), equilibrated and then subjected to a heat pulse. From the results shown in Figure 
13 we conclude that a small heat pulses does not change the state of the array at 9 <C 9q, while a large one brings 
it above the effective glass transition temperature and eliminates all memory of its initial state. If the metastable 
states did subdivide with decreasing temperature, the correlation and response functions obtained with and without 
the application of a moderate heat pulse would be different; in the former case memory of some of the subdivisions 
would be erased. Clearly this is not the case for the periodic array, and response functions for this cooling-heating 
process (not shown) are also consistent with this conclusion. Furthermore, we checked that the results obtained 
when the system was cooled slowly to 9 = 1.4, equilibrated and then subjected to the heat pulse were qualitatively 
similar. ^From these results we conclude that there is no subdivision of metastable states below the glass transition 
temperature. 
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Fig 13. The correlation function, d(t,t'), for the cooling-heating regime. Here the temperature was reduced infinitely fast to 
8 = 1.4 at t' = 0; the system was equilibrated for < t' < 30, then heated linearly for 30 < t' < 40, then cooled linearly back to 8 = 1.4 
for 40 < t' < 50 and then measured at t = 70. Different curves correspond to different amplitudes of the heat pulse as shown. We see that 
for heat pulses of amplitude less than 0.6 the system recovers its original state; however for larger heat pulses all memory is lost. 



We also checked that these results were not due to the leading-order approximations made in the derivation of the 
dynamical equations, and thus performed direct simulations of the array to confirm that this conclusion remains valid 
at low temperatures and at large a. First we simulated finite-size arrays with a = 1.0 at very low temperatures and 
confirmed that the number of metastable states grows exponentially with system-size. We also studied the distribution 
of overlaps between these states, and our results are displayed in Figure 14 for a = 1.0. Here we see that the overlap 
distribution evolves to a sharp peak with a maximum that shifts to lower q with increasing number of wires N. In 
particular we do not observe any states with large overlaps for N > 15; furthermore P(q) suggests that for large 
values of N the states are equidistant in phase space. We also found qualitatively similar results for smaller values of 
a. 
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Fig 14. Distribution of overlaps for finite-size arrays with a = 1.0. Here the gray, dashed and black curves correspond to 
N = 7, 15 and 31 respectively. 



The equidistance of metastable states does not imply that the structure of the phase space is simple. In particular 
the basins of attraction associated with different metastable states appear to be complicated. For example, we have 
studied the probability, Q(q), that a randomly chosen state that has overlap q with a metastable one evolves to this 
particular metastable state. The results are displayed in Figure 15. We see that the probability to evolve to different 
states remains finite even if the starting overlap is close to unity. Our results suggest that 1 — Q(q) oc (1 — q) for 
N — > oo. If true, this has important consequences for the distribution of barriers. First we note that the energy 
required to flip m spins cannot be larger than mTn. The probability of evolving to a different metastable state after 
the flipping of m spins is roughly j± so the probability to find some m spins which lead to this result is 0(1). This 
implies that even in the infinite-range periodic system there is a broad distribution of barriers which extend down to 
Eb ~ T . We note that this result implies that the distribution of barriers and overlaps are not related. 
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Fig 15. The probability, Q(q), to evolve back to a metastate state, A, starting with a randomly chosen state with 
overlap q with state A. The data shown here are the results of a direct numerical simulations of a finite size arrays with a = 1.0; the 
gray, dashed and black curves correspond to N = 7, 17 and 23 respectively. 



IV. THE PHYSICAL CONSEQUENCES 



What are the observable consequences of the different behavior in the response and the correlation functions for 
the three cooling regimes? If the periodic array is placed in a time-dependent field H (t) that is small compared to its 
constant part, the resulting Josephson currents generate a total magnetic moment 

M t = ii f g tt ,d tt , {H{t) - H(t')} dt' (39) 

J OC 

which can be measured with an inductance coil; here 



where L — Nl. The form of M in (39J) should be contrasted with the response of the spin systemto a physical magnetic 
field which is just g t t'\ the dissimilarity is due to the fact that in the Josephson array direct coupling of the physical 
field to the "spins", Si = e 1 ^, is prohibited by gauge invariance. Instead H couples to the network via changes of the 
Josephson exchanges; hence the response involves a four-spin term (GD). Here we discuss the magnetization induced 
by a small increment of magnetic field applied at time this is of course the simplest possible case. Still the value 
of the observed moment will depend on the relative magnitudes of tn, t Cl t* and t as indicated in the schematic of 
Figure 16. We proceed to discuss each separate case below. 
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t f = (t w + t c ) 




time 

Fig 16. Schematics of (a) Temperature (T) (b) Application of an Additional Field (Ai/) vs time indicating the time-scales 
involved in a finite-cooling experiment as discussed in the text; Tq and tQ refer to the temperature and time when the system goes out 
of equilibrium, and tu and t* are the time-scales associated with the onset of an additional field and ageing; t is the total time associated 
with a given measurement. 



A. t,t H > t c 



Here we apply a field after a fast quench, a regime where there is a generalized fluctuation dissipation relation. 



Then the expression for the moment in (39) becomes 



M = iix f" \{D 2 )'AHdt' = ^xd 2 7 &(t H ~ t G ) AH (41) 

Jt c ~o 2 2 V t-ta J 

where we have used the "fast cooling" scaling form for the correlation function exhibited in Table I. Since the fast- 



cooling response function decays with time, (41) indicates an "ageing" moment which decreases as a function of 
decreasing ratio t f~ t ta (tit it > tG). ( fll| ) also implies markedly different behavior if tjj > tQ or < tc', this is the 
array analogue of the standard zero-field vs. field cooled susceptibility measurement often performed on spin glasses.0 



B. t H < t c = t 



Now the cooling is continued after the field application (cf. Figure 13), and the dominant contribution to the 
expression for the moment, (|39|), will come from the range t c — t* < t' < t c where t* is the "ageing" time-scale 
discussed in the slow cooling regime in Section III. Here there are two possibilities: 

(a) t c -t H < t* 

In this case the correlation and the response functions are non-decaying (cf. Figs. 4-7) and the moment is a constant 
as in (||). 

(b) t c -t H > t* 
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The correlation and response functions have no memory on time-scales greater than t*; therefore the moment is 
zero. The system is displaying ageing analogous to that observed in spin glasses;0 as in the latter case, there appears 
to be a "rejuvenation" associated with improved memory at low temperatures where t* increases. 

We also remark that application of a a.c. held could lead to a direct experimental determination of G W D W ; there 
we expect ~ fi^t*) as shown below in Figure 17. Here we define = jl^, dt' exp 11 ^* - * ) dtt'Qtv ■ 

We note that the low-frequency part of the displayed curves is universal and so would provide an experimental 
determination of the time-scale associated with ageing, t* . 




Fig 17. Typical imaginary parts of ac susceptibility for the slow cooling processes. Here we show the results for the cooling 
rates T = 0.00057 (upper curve) and T = 0.0022 (lower curve). The upper curve shows imaginary part of susceptibility measured at 
reduced temperature 8f = 1.0; the lower curve corresponds to the measurement at 8f = 1.2. The frequency w is measured in units of the 
inverse characteristic time, t* = 120, corresponding to the decay of the correlation function. Here we define it by the following procedure: 
(i) we do a linear fit to d tt i ; (ii) we identify the intersection of this approximated version of d tt i and the x-axis as t — t* . Using this 
definition, t* acquired the identical value for the two cooling processes considered here. 



C. t > t c > t H 



In this "crossover regime" the incremental time-dependent field is applied before the system has reached its final 
cooling temperature (Tf < Tq)\ then the expression for the moment becomes 

M - / TG[f(t - t')]D[f(t - t')]AHdt' = ha^H (42) 
Jt G 

where [ia is a number; this result is obtained by changing the integration variable to t = ft' where f ~ (^f) 1 v< ~ Tf \ 
In this regime the moment is simply a constant with no time-dependence but with possible dependence on the cooling 
rate if 77 =^ 0. 
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V. CONCLUDING DISCUSSION 



In conclusion we have studied history-dependence and ageing in a periodic model in the absence of disorder. In 
the case of a fast quench (t ^> t c ) a generalized fluctuation dissipation theorem relation (G ~ xD') exists between 
the correlation and the response functions, and furthermore the latter exhibit a predicted scaling formJia In the other 
limit of slow cooling [t = t c ) there is numerical evidence for the existence of an "ageing" time-scale t* such that 
for t > t* the system forgets its past, but for t < t* it has "perfect memory". This additional time-scale increases 
to that associated with the cooling process with decreasing temperature, indicating that at low temperatures the 
system recovers its memory and an adiabatic solution.-(cf. Table I) reappears. However the form of this solution 
does not ressemble (cf. Fig. 8) the "adiabatic ansatz"ll3 G ~ S and D ~ Dq, which is probably never correct. We 
speculate that this is due to the presence of an extensive number of states exactly at the glass transition, where the 
time-scale associated with the cooling process is significantly less than that of the system's relaxation. By contrast, 
in the Sherrington-Kirkpatrick model such "branching of states" occurs at lower temperatures wher.e-.the cooling and 
relaxation times are comparable; in this case the assumption that there are no additional time-scalestij leads to results 
which agree with those obtained by thermodynamic methods.Ej 

We also speculate on two different scenarios for the decrease of the exponent 77 (t* = (^r) 1 * ? ' T ' ) ) with decreasing 
temperatures, and the reappearance of the adiabatic solution. One is that the latter solution is always present for 
9 > 9g, but its weight is proportional to a high power of 6 — Og, and thus is only observable at lower temperatures. 
Another possibility is that 77 — >■ corresponds to a second dynamical transition at lower temperatures. Below the 
temperature associated with this second dynamical instability the metastable states would not evolve dramatically 
with subsequent cooling, leading to the absence of further ageing effects. Naturally the emergence of t* from an analytic 
ansatz for the dynamical equations would provide additional insight to this issue. Unfortunately these equations are 
difficult to solve in the slow-cooling regime since the corrections to long-time behavior are of the same order as the 
leading terms themselves. 

Just below its glass transition temperature, the periodic array becomes stuck in one of an extensive number of 
metastable states that are equidistant in phase space. Our numerical studies indicate that there is no further subdi- 
vision of these states at lower temperatures, in contrast to the situation for the (disordered) Sherrington-Kirkpatrick 
model. However the probability to evolve to different metastable states remains finite even if the starting overlap is 
very close to unity, which suggests that the distribution of barriers and the overlaps are not simply related. It would 
be interesting to know if these features of the periodic array's low-temperature phase space are signatory of a glass 
with distinct dynamic and static transitions. 

Naturally there remain several questions to be addressed in future projects. A key feature of this periodic array 
is that it can be realized experimentally;!!] a calculation and subsequent experimental determination of its nonlinear 
response be quite interesting. In a real system there arc evcrpresent induced screening currents ignored in the present 
treatment, and the addition of these terms could lead to novel memory effects. Amusingly enough, ageing effects are 
not present in the disordered long-range array; a study of the crossover between the periodic and the random networks 
might provide further insight into the origin of this widespread but mysterious phenomenon. 
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